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Abstract. A discontinuous Galerkin method for the ideal 5 moment two-fluid plasma 
system is presented. The method uses a second or third order discontinuous Galerkin 
spatial discretization and a third order TVD Runge-Kutta time stepping scheme. The 
method is benchmarked against an analytic solution of a dispersive electron acoustic 
square pulse as well as the two-fluid electromagnetic shock [1 1 and existing numerical 
solutions to the GEM challenge magnetic reconnection problem |2|. The algorithm can 
be generalized to arbitrary geometries and three dimensions. An approach to marn- 
tarnrng small gauge errors based on error propagation is suggested. 
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1 Introduction 

Fusion power promises to be a safe, efficient and environmentally friendly energy source. 
Controlled fusion power concepts have been under investigation for decades, the vast 
majority of these concepts require an intimate understanding of plasma physics to deter- 
mine the stability and confinement properties. Numerical plasma physics has proved ex- 
tremely valuable in deciphering experimental data and predicting the behavior of plasma 
experiments. Many plasma fluid models, and in particular the full two-fluid plasma 
model, have received very little attention from the numerical plasma physics community. 
This work describes an advanced algorithm for the ideal 5-moment two-fluid plasma sys- 
tem. 

To solve problems in plasma physics and to gain physical intuition of plasma phe- 
nomena a hierarchy of classical plasma models have been developed. The most funda- 
mental continuum plasma model is the Vlasov model which eliminates individual par- 
ticles in favor of a continuous distribution function. This model is six dimensional as 
the distribution function is a function of both position and velocity. The Vlasov model 
can be re-written as an equivalent system that consists of an infinite number of moment 
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equations. A reduction of the Vlasov model can then be obtained by truncating this in- 
finite series. Assuming scalar pressure and setting the heat tensor and higher moments 
to zero produces the 5 moment truncation of the Vlasov model. This model is known as 
the ideal 5 moment two-fluid plasma model, and will be discussed in this paper. Asymp- 
totic approximations of this two-fluid system produce a series of increasingly simpler 
fluid models including two-fluid MHD (Magnetohydrodynamics), Hall MHD and then 
the ideal MHD models. 

The main benefit of a fluid model over the Vlasov model is the reduced dimensional- 
ity from 6 dimensions to 3 dimensions. Physics is lost in this reduction, but an enormous 
amount of physics relevant to fusion and spacecraft propulsion remains in the fluid de- 
scription. Ideal MHD has been extremely successful in explaining large scale instabilities 
in such devices as the Z-pinch, spheromak and tokamak [3.4|. Unfortunately there are 
many regimes where the description is invalid and where it fails to explain the observed 
phenomena. An example of this includes ion demagnetization which is important in 
Field Reversed Configurations [5| and Hall thrusters. Hall MHD addresses both these 
issues but fails to describe other plasma phenomena such as the demagnetization of elec- 
trons in regions of low magnetic field which is important in collisionless reconnection. 
The two-fluid MHD approach adds terms such as electron inertia which is an important 
mechanism for breaking the frozen in flux condition for electrons as it acts as a "dissi- 
pation" mechanism in the absence of resistivity ||6|. The quasi-neutrality condition still 
constrains the electron and ion motions, to allow complete independence of electron and 
ion motion the quasi-neutrality condition must be relaxed; the result is the ideal two-fluid 
plasma system. 

Two-fluid effects are important in the generation of turbulence through microinstabil- 
ities. Most plasmas are turbulent at some scale, however the simplest fluid model, ideal 
MHD, describes plasmas physics that is more or less laminar where the two-fluid model 
produces turbulent phenomena. This can be explained in part by the fluid description 
of electrons. In a two-fluid model both the electrons and the ions may become unstable 
independently. In particular, electrons carry most of the current in an MHD plasma. This 
current may produce a large amount of differential motion in the electron fluid when 
magnetic field gradients are present even if the plasma is in a static MHD equilibrium. 
The generation of microturbulence through processes such as the lower hybrid drift in- 
stability and the modified two-stream instability may be important in both Z-pinch and 
theta-pinch plasmas. These instabilities are frequently cited as sources of anomalous re- 
sistivity 1 7 1, magnetic diffusion and heating and in certain cases may ultimately drive 
macroscopic MHD instabilities [8|. 

A particularly good application of the two-fluid plasma model is the fusion Z-pinch 
Q. Many plasma experiments last a few seconds whereas the shortest plasma times 
scale, the electron plasma oscillation, can occur on the scale of pico seconds. However, in 
the case of the fusion Z-pinch these time scales can be compressed to about 4 orders of 
magnitude between the shortest time scale, the electron plasma period the MHD insta- 
bility growth time, which puts two-fluid Z-pinch simulations in the range of numerical 
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methods. Conceptual Z-pinch fusion reactors are high density with extremely strong 
magnetic fields which makes the two-fluid plasma system particularly applicable fW\. 
Furthermore the radius of the pinch can be 1000 Debye lengths in some designs 
which means Debye length scales could be resolved in very high resolution simulations. 
Artificially increasing the electron mass to ion mass ratio and increasing the ratio of the 
Alfven speed to the speed of light can make the Z-pinch problem more computationally 
tractable while maintaining the relevant physics. The analysis of microinstabilities such 
as the lower hybrid drift instability may be important in progressing towards a better 
understanding of Z-pinch physics. 

Algorithms have been designed for various fluid plasma models including MHD 
p2p3 L Hall MHD 1 14^171 various forms of electrostatic two-fluid plasma models 118^20) 
and the ideal two-fluid syste m |[l 21 22) . The first well described two-fluid plasma al- 
gorithm was the ANTHEM |j21||22| code. It was used to simulate fast phenomena in 
high density plasmas inaccessible to PIC codes, the applications included simulations of 
plasma opening switches and the fast igniter concept. ANTHEM used a flux corrected 
transport (FCT) type algorithm for the fluids and an FDTD type algorithm for the fields. 
This type of algorithm is difficult to extend to general geometries because of the staggered 
scheme used for Maxwell's equations. In 1 1 1 a full two-fluid algorithm using the finite 
volume method for both fluids and fields was described for one-dimensional problems. 
An issue with the finite volume algorithm was the decay of equilibrium solutions due to 
the low order of accuracy which resulted from the source term integration and diffusive 
limiters used. Major improvements on the finite volume technique have been made with 
the help of various divergence cleaning techniques and careful attention to source term 
treatment. An improved finite volume approach is described in [231. 

The purpose of the paper is to develop a numerical algorithm for the ideal 5 moment 
two-fluid plasma system using the discontinuous Galerkin method so that it can be easily 
generalized to arbitrary geometries and to arbitrarily high order accuracy to help capture 
plasma instabilities. TVB discontinuous Galerkin methods are described in [24-26]. They 
are extended to the multi-dimensional Euler equations in p7| and to Maxwell's equations 
in 1 28 29\. The two-fluid system of equations consists of two sets of Euler equations, one 
for the electrons and the other for the ions, and the complete Maxwell's equations. The 
ideal two-fluid system differs from the ideal MHD equations in that it is composed of 
three separate (but well understood) hyperbolic systems coupled through source terms. 
The MHD equations are, on the other hand, a unique hyperbolic system. A discontinuous 
Galerkin method for the MHD equations was developed in |30| and for two temperature 
MHD in |31|. The technique is used in a Vlasov-Maxwell algorithm in [ ,32] , numerous 
other applications can be found in p3| . 

In section |2] the ideal 5-moment two-fluid model is described and the equations are 
presented. In section |3] a scalar model problem is derived from the two-fluid systems 
which helps to illustrate some of the numerical issues with the system. In section |4] the 
discontinuous Galerkin method applied to the two-fluid plasma system is presented. In 
section |5] simulations are presented including an electron acoustic pulse for validation 
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of code accuracy, an electrostatic shock, the two-fluid electromagnetic shock |[TJ, and the 
GEM challenge magnetic reconnection problem [2| where the reconnected magnetic flux 
can be compared to published results. Finally, in section |6]the conclusions are discussed. 

2 Two-Fluid Model 

The full two-fluid plasma model consists of a set of fluid equations for the electrons and 
ions plus the complete Maxwell's equations including displacement current. The fluid 
and electromagnetic systems are coupled by Lorentz forces and current sources. In the 
following equations, E is the electric field, B is the magnetic field, qs is the species charge 
(subscript s is / for ions and e for electrons), ps is the species density, nis is the species 
mass. Us is the species velocity, Ps is the species pressure and eg is the species total energy 
with Cs = jPsUg + jt^ps- The species number density is defined as Ug = eo is the 
permittivity and f/o is the permeability of free space. Maxwell's equations are presented 
in SI units. The complete Ampere's law is used 

dtE-c\VxB) = --J2—PsUs, (2.1) 



and the complete Faraday's law 

The magnetic flux equation, 
and Poisson's equation ( |2.4| |, 



9fB + (VxE)=0. (2.2) 



V-B = (2.3) 



V-E = -{qini+q,n,) (2.4) 
^0 



are constraint equations which can be derived from Ampere's law ( |2.1| , Faraday's law ( |2.2| 
as well as the species continuity equation given below in ( 2.11| under the assumption that 



the constraints are satisfied initially. A simple way to reduce the error in the constraint 
equations is to use the perfectly hyperbolic Maxwell's equations where equations ( |2.1| |, 



|2.2| |, (2.3 1, ( |2.4| | are modified, so that Ampere's law becomes 



dtE-c^VxB) + VxpE = --T^PsUs, (2.5) 

Faraday's law becomes 

afB + (VxE) + Vi/^B=0. (2.6) 
The magnetic flux equation becomes, 

1 dxpE 



r| dt 



+ V-B = (2.7) 
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and Poisson's equation becomes. 



1 ^'^^+V-E = -{qini+qene). (2.8) 



Where tpB and tpE are correction potentials and Fg, Te are the correction potential wave 
speeds which propagate errors in the divergence constraints out at the speeds Fg and Te- 
The correction potential wave speeds can be set to zero for problems where the correction 
terms are not necessary. Details of the correction potential technique are described in 



1 34 1 and used in our previous work p3| . In many problems of experimental interest 
electrons flow from a surface into the simulation domain. In these situations we believe 
more sophisticate divergence preservation will be required. It's important to maintain 
exact charge conservation especially when the emission is space charge limited. For this 
type of problem a constrained transport approach should be taken. Charge conserving 
constrained transport is frequently used in particle in cell (PIC) codes |35| , [ [36| , [371. 
In addition, V • B = preserving constrained transport is frequently used in the MHD 
system [38j, (39), ||40|, ||4T). 

The fluid equations are simply the inviscid Navier Stokes equations with Lorentz 
force source terms. Each fluid species has its own equation for energy, 

dtCs + V- {Us {es+Ps)) = — psE • U„ (2.9) 

momentum, 

dt {pslJs) + Vcc (psU^U,) + V« (s^f^Ps) = ^ps (E+U, X B) , (2.10) 

\ / tils 

and continuity, 

dtps + V-{pslJs)=0. (2.11) 

This means that each species has its own temperature, velocity and number density. As 
a result, quasi-neutrality is not assumed and things like electron plasma waves and ion 
subshocks should be observed numerically. This system is identical to the system used 
iny. 

The ideal two fluid plasma system can be written as three systems of balance laws. 



for the electron equations. 



dt 



dQi 



+ y-Fe{Qe)=^e{Qe,Qem), (2.12) 



dt 



^y-Fi{Qi)=Xpi{QuQem), (2.13) 



for the ion equations, and 



dt 



+ y-Fem{Qcm)=tpem{Qi,Qe), (2.14) 
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for Maxwell's equations. These balance laws, Eqns.( Z12| -( 2.14[ |, are given in full form by. 
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3 Derivation of a Scalar Model Problem 



The ideal 5-moment two-fluid system is unusual in that the source terms act as harmonic 
oscillators, the source terms are purely dispersive without dissipation or amplification. 
This fact is important when considering numerical methods to use, as the source term 
integration should introduce as little dissipation as feasible in order to avoid damping 
these oscillations. It frequently occurs in a two-fluid plasma that convective forces are in 
balance with oscillating sources to produce an equilibrium. With this in mind, a simple 
model problem is derived which may help one choose a proper numerical method for 
the two-fluid system. 

Linearizing the electron x-momentum equation ( 2. 10| and Ampere's law ( |2.1| | while 
assuming a constant background ion density and assuming that B = 0, a partial differen- 
tial equation for electron plasma oscillations can be derived which takes the form 

J^ = -^peU, (3.1) 
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where u is the perturbed x velocity and cOpe is the electron plasma frequency. This can be 
transformed to a first order equation in complex variables by making the transformation 
ajpeV = ^ and letting Q = v + iu. 

^ = icOpeQ. (3.2) 



By transforming variables Q{x,t) — >• Q{r],t) where rj = x-\-at, Equation( |3.2| | becomes an 
advection oscillation equation 



dt dt] 



+a^ = icOpeQ, (3.3) 



where a is the wave propagation speed and cOpe is the oscillation frequency. Any algo- 
rithm that is stable to the two-fluid system should be stable when applied to Eqn( |3.3| |. 
This equation has solutions Q = Ae'^*^''^'^'' where (v = ak—cOpe. In particular( |3.3| l admits a 
steady state solution Q = Ae''"' on an infinite domain where the source term is in balance 
with the flux. It is important to note that there are an infinite number of equilibria which 
differ by a continuous range of scalar factors A. This point is important when considering 
numerical methods for this system since a numerical method with too much dissipation 
could conceivably move a steady state solution from one equilibrium to another all the 
while moving towards a state where A = 0. This loss of amplitude is also observed in 
equilibrium type problems in the two-fluid system. An effective numerical algorithm 
must be both stable to the advection equation and oscillation equation and must have 
low dissipation in equilibrium type problems. The discontinuous Galerkin method is an 
ideal candidate for solving the two-fluid system. Its accuracy can easily be increased to 
reduce numerical dissipation while being stable to both the advection equation and the 
oscillation equation. 

4 Runge-Kutta Discontinuous Galerkin Method 

Discontinuous Galerkin methods are high order extensions of upwind schemes using a 
finite element formulation where the elements are discontinuous at cell interfaces. Details 



of the method are discussed in |24-27| and reproduced here for our particular case. The 
balance law 

^ + V-FiQ) = ip{Q) (4.1) 

is multiplied by the set of basis functions {vr} and integrated over the finite volume 
element K. For second order spatial accuracy the basis set on a unit square reference 
element is 

{vr} = {vo,Vx,Vy} = {l,x,y}. (4.2) 
For third order spatial accuracy 

1 1 

{vr} = {vo,Vx,Vy,Vy:y,Vxx,Vyy} = {l,x,y,xy,x^ - -,y^ --} (4.3) 
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is used. The equation is written. 



[ -^VrdV+ [ {V-F)v,dV= [ xpVrdV. (4.4) 
Jk ot Jk Jk 



Integrate by parts to get 



I -^VrdV+ I (F-n)vrdT- [ F-(Vvr)dV= [ itv,-dV. (4.5) 
Jk at JdK Jk Jk 

The surface flux f ■ n is a numerical approximation of the exact flux F ■ n across the inter- 
face. The discrete conserved variable Q is defined as a linear combination of the basis 
functions inside an element K, with 

Q = Y2vrQr. (4.6) 

r 

The integral J^^VrdV = ^CV where C is the constant J^v^dV and V is the volume 
of the element. Using these definitions we get the discrete equation 

^Cy + ^^W; [Fi-n] VrlTe-J2^mFm - {^Vrm)V = J2^m^mVrmV , (4.7) 
"''el m m 

when the integrals are replaced by appropriate Gaussian quadratures. is the surface 
area of the cell face in consideration, e refers to an element face, / are quadrature points 
on a face with Wi the associated weight, m refer to quadrature points in the volume with 
Wm the associated weight. Functions with subscript / or m are evaluated at the face 
quadrature points and m^^ volume quadrature points respectively. For a second order 
method the edge integrals are replaced by a two point Gaussian quadrature 

/_;/W.../(-L)+^(__L) ,4.8, 

A four point quadrature is used for the volume integral given by. 



/j'/(x,y)«y./(-L,-L)+/(--L,i.) 



The discrete equations for the second order scheme using the basis functions given in ( |4.2| 
become 

^y + ^^W/ {Fi-ne)voiTe = J2^'n^mVOmV, (4.10a) 
"'el m 

^^y + 3^^ro/ {Fi-ne)v^iTe-3j^WmFm-{VVxm)V = 3j^WmXpmVxmV, (4.10b) 

e I m m 
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3Qi 

dt 



T e-'^Y^WmFm - {V Vym) V = 3j^WmrpmVyntV ■ (4.10c) 



e I 



The derivatives of the basis functions can be calculated analytically since the polynomial 
basis functions are known. The discontinuous Galerkin method is applied to each bal- 
ance law ( 2.12"| |(2.13l( 2.14| | at every time step. For the third order space method the edge 
integrals are done using a 3 point quadrature 



V3\ r( 



(4.11) 



The volume integrals are performed using a 9 point quadrature which can be calculated 
by doing a 3 point integration in the x direction and then a 3 point integration in the y 
direction. This produces the following approximate integral 



25 
81 



/(^,f)./(-^,f)./(-^,-f)./(^,-f 

/H)^/(#o)./(o,-f)./(-^ 



+ 



40 
81 



,0 



(4.12) 



For the 3rd order scheme the following discrete equations must be updated in addition 



to those given by the second order scheme (4.10 1 using the basis functions defined in ( |4.3| 
becomes. 



^^y + 9^^W/ {Fi-He) V^yll e-^Yj^mFm - {'^ Vxym) ^ = 9^ l/^m^^xy/ ^/ (4.13a) 



e I 



^Q^xy j^'^YJ£^Wi {Frne)Vx^iTe-^Y^WmFm-{'^V^xi)V=^Y2wmtpmVxxmV, (4.13b) 



3^ 



e I 



{Prne)VyylTe-^J2^mFm-{yVyym) V =^Yl^m^'nVyymV . (4.13c) 

Though the spatial discretization uses a finite element approach, the time integration 
uses standard finite difference methods which are described in the next section. The 
algorithm described is an explicit finite element method, data is only exchanged between 
neighboring cells. The solution does not need to be continuous at cell interfaces which is 
particularly useful for problems with shocks. 
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4.1 Time Integration Schemes 

Time integration schemes that are stable for the advection equation must also be stable to 
the oscillation equation if they are to be stable in general to the two-fluid system. In this 
paper the 3rd order TVD Runge-Kutta method |24J is used. 



Qi = Q"+AfL[Q"] 



(4.14a) 



Q 



lQ" + l{Q'+AtL[Q^]), 



(4.14b) 



(4.14c) 



The time integration scheme is applied to each at every time step to evolve the solu- 
tion. The term L[Q"] represents the entire "left hand side" which is everything but the 
time derivative evaluated at Q". It is important to note that all two step second order 
Runge-Kutta schemes are unstable to the oscillation equation p2|. 



4.2 Evaluating F • n 

The flux F ■ n can be evaluated a number of different ways. The local Lax flux is used in 
this paper and is computed at each face as 

F-n = l{Fi'+Fr^,)-n-^-\Mi^^/2{Qt-Q;_,^)-n, (4.15) 

where | A | /2 is the maximum eigenvalue of the particular system based on the averages, 
Qo, of the conserved variables at the centers of cell / and i + The local Lax flux is a well 
known flux function that can be used in the discontinuous Galerkin method ||24|. For the 

fluid systems | A | ,+1/2 = ( | | + ( 7a ^ ) ^ ) is used. For Maxwell's equations | A | = c. 

The superscripts + and — mean that the Q is evaluated at the upper or lower edge of the 
cell. 



4.3 Limiting 

High resolution schemes typically use limiting to prevent spurious oscillations near dis- 
continuities and for stabilization of non-linear systems [43 1. Limiters can also be used in 
the discontinuous Galerkin method, though instead of being TVD, minmod limiters pro- 
duce a scheme that is TVDM or TVD in the mean. This means that the solution is TVD in 
Qo, but not necessarily in Q. 

Following the procedure described in f25\ the conserved variables Q can be limited in 
terms of characteristics or in terms of components. To limit Q in terms of characteristics 
the Q are first transformed to characteristic variables g where g = LQ and L is the left 
eigenvector matrix of the flux Jacobian calculated from Qq. The left eigenvector matrix 
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is also applied to the differences L [Qq - QqJ = A+^o and L [Q'^ - Qq ^ j = A^^o- Limit- 
ing is performed directly on transformed variables and then the solution is immediately 
transformed back to determine the limited form of Qx, 



Qx = L ^m{gxA^goA go) 
where m is the minmod limiter defined by 



(4.16) 



m{a,b,c) 



' max(a,b,c) if sign(fl) =sign(&) =sign(c) : 
mm[a,b,c) if sign(fl) =sign(&) =sign(c) : 
otherwise 



(4.17) 



The minmod limiter Eqn. ( 4.17| | is typically used for each of the fluid equations Eqn.( 2.12| ( 2.13| l 
while the modified minmod limiter can be used to reduce the dissipation 



m{a,b,c) 



'a \{\a\<Mdx'^ 
m[a,b,c) otherwise 



(4.18) 



where M is a constant. Component limiting is done in a similar manner except no 
transformation is necessary, so that the limiter is directly applied to the variables Q. 
Component limiting has the advantage that it is faster than characteristic limiting and 
it does not introduce machine precision errors that can result during the transforma- 
tion Q = L^^ (i-Q)- The disadvantage of component limiting is that the approach is not 
TVDM |25| and unphysical oscillations can appear in the solution. In this paper charac- 
teristic limiting is used. 

When a 3rd order DC method is used, two types of limiters can be used. The first 
method follows the procedure of second order method and is described in |26| where 
if Qx 7^ Qx then all higher order coefficients are set to zero. This method is simple to 
implement and is used in this paper. 

A different and potentially better 3rd order limiter is that of |44|. In this method, the 
linear terms, Qx, Qy, Qz are limited in the same way as the second order method while 
the higher order terms Qxy, Qxx and Qyy are limited as follows. 



Q 



1 

XX ' 



and 



Q 



-m 



■ m 



Qxx' 2 ( Qx 



i+l 



-Q'x 



Q 



-Qx-' 



Qyy'' 



(4.19) 



(4.20) 



finally, the term Qxy is limited by setting it to zero if either Qxx or Qyy is limited. In f26l 
it was suggested that the high order terms could be limited by simply setting them to 
zero if the linear terms are limited. The justification is that oscillations in the higher order 
terms would only be important when oscillations in the linear terms exist. 
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4.4 Stability 

The stability limits of the numerical algorithm just described are defined by the highest 
oscillation frequency of the system or by the CFL condition based on the speed of light. 



Typically the highest oscillation frequency is the electron plasma frequency Wpe= 



and a time step is chosen for which the time integration scheme is stable to this frequency 
of oscillation, this time step is typically At<^. When the CFL condition dominates the 

time step < g ^ is used for the second order spatial discretization in 2D and < ^ ^ 
for the third order spatial discretization in 2D. 

5 Simulations 

The two-fluid system describes many dispersive waves including the electron acoustic 
wave. In a plasma the electron acoustic wave is coupled to the plasma frequency produc- 
ing a wave that is essentially stationary for sufficiently long wavelengths or sufficiently 
cold plasmas. In the following a dispersion relation for the electron acoustic wave in a 
warm plasma is derived from linearized two-fluid equations. The dispersion relation is 
used to calculate an analytic solution to the propagation of an approximate square pulse 
in a two-fluid plasma. The numerical two-fluid solution using the 2nd and 3rd order 
discontinuous Galerkin method is compared to the analytic solution using various grid 
resolution. The order accuracy of the algorithm in the L2 norm is calculated from these 
results. 

In real kinetic plasmas damping such as Landau damping may result in the damping 

of waves observed in simulations such as this. These results are meant to illustrate veri- 
fication that we are solving the correct system in addition to showing we can achieve the 
desired accuracy. 

Assume infinitely massive ions with a backgroimd number density no for both elec- 
trons and ions and charge qi = —q^- Furthermore, assume background electron and ion 
pressures Pq while all other background quantities are zero. A perturbed electron ve- 
locity u] = Uoe'^^"^^^"*^ is assumed. Corresponding perturbed electric field, density and 
pressure profiles can be derived from Poisson's equation, the continuity equation and the 
energy equation so that the perturbed electric field = ^^noqeu], perturbed electron 

pressure, Pg„ = ~ (^^^JePoul, and perturbed electron density, pl„ = — (^^^ Poeu]- The 
electron acoustic dispersion relation is 





(5.1) 



The positive root defines waves that travel to the left. A square pulse on a periodic do- 
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main is defined by taking linear combinations of these waves, 

CO ■ 

n=0 

As a practical matter, an approximate square pulse is used since the high wave numbers 
cannot be resolved numerically unless the spatial resolution is sufficiently high. Figures 
[T]and|2]illustrate this dramatically. Figure [Tjshows the initial square pulse initialized with 
5000 wave modes. Figure [Tjshows the analytic solution after t = 1000 steps. Since there 
is no dissipation in the system and the system is dispersive the high frequency modes in 
the initial conditions play an important role in the final solution; this makes the issue of 
convergence in shock-type problems that start out with discontinuities difficult to assess. 
As a result, in the simulations and analytic solutions that follow, n = 9 is the highest 
mode included in the expansion and k„=2nn. Finally, only the real part of all perturbed 
quantities are used in the initial conditions, thus, 

El{x,t) = Y,El„ (5.3) 

n=0 

P}ix,t) = f^P}„ (5.4) 

n=0 
9 

n=0 

In this simulation qi = —qe = 10, Cq = 1, = 1, Pq = 1, = 1, m, =oo, 7^ =2 for convenience. 
To ensure that the solution is in the linear regime Uq must be set to a small value. For 
these simulations Uq = 1x 10^^. The domain is periodic with length 1. Electromagnetic 
waves do not exist in this problem, but the speed of light c = 1. The simulations are 
run to time t = 3 at several different resolutions and 20,000 time steps are taken for the 
highest resolution simulation which has 320 cells. The initial conditions for the electron x- 
velocity perturbation are shown in figure|3j An approximate square wave is used to excite 
several wave modes to test the algorithms performance effectively. A problem in the 
linear regime is used for two reasons. First of all, analytic solutions exist and are easy to 
calculate, secondly, in the linear regime the limiters can be turned off so the solution can 
be observed without the added dissipation which can reduce overall accuracy making it 
much easier to compute numerically the accuracy of the algorithm. Ultimately, problems 
in the non-linear regime, problems with large Uq for example, require limiters and the 
overall accuracy of the numerical solution is reduced. The design of effective limiters 
may be the most important problem in gaining computational efficiency from 3rd order 
or higher discontinuous Galerkin methods for the two-fluid system when compared to 
the 2nd order method. 
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ure 1: Electron velocity of a square electron acoustic square pulse at time t 
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Figure 2: Analytic solution of the electron velocity of a square electron acoustic pulse at time f = 1000. This plot 
illustrates the dispersive, non-diffusive nature of the two-fluid system that makes it numerically difficult even in 
the linear regime. The high frequency modes still contribute significantly to the amplitude of the solution late 
in time and they always will because there is no physical diffusion. The dispersive, non-diffusive nature of the 
two-fluid system can make it appear that our numerical solutions suffer from significant numerical dispersion 
errors when they are actually correctly capturing dispersion described by the model. 
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After 3 time units the square pulse shape has disappeared due to wave dispersion. 
Plots of the numerical solution verses the analytic solution at various grid resolutions are 
shown in figures |4| |5| and|6]at 3 time units. In figure |4] there are 40 cells in the domain 
and the 3rd order method shows evidence of resolving the highest order mode. The 2nd 
order method only captures the bulk features. In figure |5] there are 80 cells in the domain 
and the 3rd order method captures the amplitude of the highest modes and matches the 
analytic solution very well. The 2nd order method still struggles to resolve the highest 
mode (note the solution at the two skinniest spikes). In figure |6] there are 160 cells in the 
domain and the 2nd order method still does not match the amplitude of the highest order 
mode and does not match the amplitudes much better than the 3rd order method at 40 
grid cells. Figure [7| shows a plot of the convergence history of the numerical solutions 
along with calculated order of accuracy. In both cases the calculated order of accuracy 
varies from less than 1 for low resolution where the high frequency modes are barely 
resolved to better than the order of the scheme when the solution is nearly converged. 

The 3rd order method performs substantially better than the 2nd order in this lin- 
earized problem when limiter are not needed. In particular, the 3rd order method pre- 
serves amplitude much better than the 2nd order, this same phenomena has been ob- 
served in certain equilibrium type non-linear problems. 
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Figure 3: Exact electron x velocity at t=0. Initial conditions are chosen so that all waves travel to the left as 
time increases. Since the waves are dispersive the initial "square wave" disappears when the solution is allowed 
to evolve. 
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Figure 4: Numerical solution using 2nd and 3rd order discontinuous Galerkin spatial discretizations with a 3rd 
order time discretization compared to the exact solution at t=3. The grid has 40 cells so the highest wave 
number mode is barely resolved with the 2nd order method. The 2nd and 3rd order solutions differ substantially 
from the analytic solution. 
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Figure 5: Numerical solution using 2nd and 3rd order discontinuous Galerkin spatial discretizations with a 3rd 
order time discretization compared to the exact solution at t=3. The grid has 80 cells across the domain. The 
2nd order solution differs substantially from the analytic solution, only capturing the lower order modes. The 
3rd order solution matches the analytic solution well. 
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Figure 6: Numerical solution using 2nd and 3rd order discontinuous Galerkin spatial discretizations with a 3rd 
order time discretization compared to the exact solution at t=3. The grid has 160 cells across the domain. The 
2nd order solution resolves the high order modes at this resolution. At higher resolution the numerical solutions 
are visually indistinguishable from the analytic solution. 
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Figure 7: Natural log of the L2 norm verses natural log of grid spacing for the numerical solutions to the 
electron acoustic wave dispersion problem. The solutions were calculated using 2nd and 3rd order discontinuous 
Galerkin spatial discretizations with a 3rd order Runge-Kutta time discretization. The numbers near each line 
give the slope of the line and hence the measured order of accuracy of the scheme. The 3rd order method gives 
substantially better accuracy than the second order method. Grid resolutions used to construct this plot are 
1/20, 1/40, 1/80, 1/160, and 1/320. The measured accuracy was computed at i = 3. 
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5.1 Two-Fluid Electromagnetic Plasma Shock 



The two-fluid electromagnetic plasma shock is an extension of the Brio and Wu shock | [45| 
to the two-fluid plasma model. The simulation was first performed in |[l}[46| and is used 
in the current paper as a benchmark. The ideal two-fluid system has no dissipative terms 
however an artificial viscosity exists due to the numerical discretization. Wave steep- 
ening of the two-fluid solution is limited by physical wave dispersion, when wave dis- 
persion is not sufficient to limit the steepening, artificial viscosity limits the steepening, 
hi real collisionless shocks in experiments and in space, kinetic effects limit the wave 
steepening when dispersion is not sufficient to limit the steepening. This simulation is 
meaningful since it illustrates the range of physics that the two-fluid system describes. 

In this paper the shock will be presented differently than in The initial discontinu- 
ity is allowed to evolve in time until the shock structure spans lOOOr^, where r^i is the ion 
Larmor radius. Time is measured in terms of light transit times across the entire domain, 
Tc = Snapshots of the shock earlier in time correspond to larger characteristic ion 

Larmor radius, ^ < 1000 where L is the span of the shock, and so the solution evolves 
from a "gas dynamic" regime of short time scales and large characteristic ion Larmor 
radius 3> 1 to an "MHD" regime of long time scales and small characteristic Larmor 
radii ^ 1. Parameters used are, qi = 10, qe = —10, and eQ = \, Ho = \, c = 1, 7^ = 7, = |, 
mi = l, me = The initial conditions on the left half of the domain are given by. 
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(5.6) 



The spatial units of figures [Sj [9} [Toj [TTj [l2j [13| are measured in ion Larmor radii r^, based 
on the initial conditions in the left half of the domain. The Debye length is \a = ^ Vgi 
also based on the initial conditions on the left half of the domain. The results in figure 
|8] correspond to a less accurate solution to that published in |[Tj figure 8. The domain of 
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figure |8]has only 500 cells in the domain corresponding to 1000 degrees of freedom (5,000 
time steps to reach this point in the simulation) while the published solution using a 
finite volume method has 4000 cells or 4000 degrees of freedom. In figure 10 the solution 
is higher resolution than the corresponding solution published in figure 9 of [1 1 as can be 
seen by the resolution of the oscillations to the left of the rarefaction wave. At this scale 
there are 5000 cells in the domain corresponding to 10,000 degrees of freedom (50,000 
time steps to reach this point in the simulation) while there are 4000 degrees of freedom 
in the previously published solution. At the final time figure 12 the solution has moved 
beyond those published and significant oscillations to the left of the shock are observed. 
Most of these waves are resolved in several hundred grid cells, and the entire domain of 
figure 12 is 50,000 cells (500,000 time steps to reach this point in the simulation). 

Due to the dispersive, non-diffusive nature of the system, the shocks solutions show 
more complexity as resolution is increased in time and space. To illustrate this figures [9| 
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13 show solutions using a time step based At = l/cOpe compared to the solution with 



At = l/ (IcOpe). These results show the large amount of damping in the At = l/cOpe case 
that results from temporally unresolved simulations, however, many of the essential fea- 
tures of the solution still remain. Similarly, increasing the spatial resolution adds higher 
frequency dispersive waves without significantly changing the overall shape. The final 
solution 12 took roughly two days on 64 processors so higher resolution solutions were 
not attempted. 
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Figure 8: Electromagnetic shocl< solution using the two-fluid equations and the MHD equations at t = O.OlTc. 
At this time the domain spans 10 ion Larmor radii or 1000 Debye lengths. It is in this regime that the two-fluid 
solution differs most significantly from the MHD or the "gas dynamic" solution. This regime has practical 
applications to Z-Pinches and PRC's due to the weak magnetization of the ions. 
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Figure 9: Comparison of 10 ion Larmor radius solution when the plasma frequency is just resolved compared to 
a solution with half the time step. Time step close to the stability limit significantly damps the high frequency 
waves virtually eliminating the trailing wave train. 
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Figure 10: Electromagnetic shocl< solution using the two-fluid equations and the MHD equations at t = 0.lTc. 
At this time the domain spans 100 ion Larmor radii or 10000 Debye lengths. Major differences from previously 
published results (Figure 9 in |1| include the large oscillations to the left of the rarefaction wave which are 
dispersive plasma waves. The differences are a result of the higher grid resolution and better accuracy of the 
algorithm used in this paper. 
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Figure 11: Comparison of 100 ion Larmor radius solution when the plasma frequency is just resolved compared 
to a solution with half the time step. 
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Figure 12: Electromagnetic shock solution using the two-fluid equations and the MHD equations at t = \Tc. At 
this point in the simulation the solution is very MHD like. However, key differences remain as a result of the 
fact that the two-fluid system models dispersive MHD waves. The most major differences are the post shock 
and post rarefaction wave oscillations. Both look numerical, but these waves are resolved in several hundred 
grid cells and are a result of dispersive plasma waves described in the model. Moving from lower resolution to 
higher resolution runs show more high frequency waves with less diffusion. This simulations was the highest 
resolution performed using 50,000 cells. In this regime the Hall MHD model should be used since if one uses 
the full two-fluid model one needs to approximately resolve the Debye length scales to get good shock solutions. 
The desire to resolve the Debye length in this simulation means that extremely high resolution was used. 
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Figure 13: Comparison of 1000 ion Larmor radius solution when the plasma frequency is just resolved compared 
to a solution with half the time step. It's very important to resolve he plasma frequency in order to get the full 
dynamics of this simulation. 
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5.2 Magnetic Reconnection 

In ideal MHD the fluid is frozen to the magnetic field lines and this prevents one field line 
from connecting with another. The addition of non-ideal terms such as resistivity results 
in fluid moving across field lines allowing magnetic reconnection to occur. Classical re- 
sistivity leads to much slower reconnection rates than are observed in collisionless space 
and fusion plasmas. Non-ideal, collisionless terms such as electron inertia and the elec- 
tron pressure gradient are partially responsible for the fast reconnection that is observed 
to occur in the earths magnetotail and fusion plasmas. In this section we show that the 
two-fluid algorithm developed in this paper produces magnetic reconnection rates that 
agree with those described of the GEM challenge magnetic reconnection problem as de- 
scribed in |2|. 

The GEM challenge magnetic reconnection problem is non-dimensionalized as in Q 



where lengths are normalized by the ion inertial length d=c/ Wpj=c yj^j time is non- 
dimensionalized by the ion-cyclotron time ^ where Bq is the magnetic field at infinity. 

The velocities are normalized by the Alfven velocity Va = (^j^^^^-f^^ ^ ■ Finally current den- 

B ZD 

sity is non-dimensionalized by Jo = and E by Eq = Bo- The domain is {—6Ad,6Ad) 
and the simulation is run out to 40/ Wd- Conducting walls are used on the y boundaries 
and periodic boundaries are used on the x boundaries. A = 0.5d, the ion to electron mass 
ratio is taken to be 25 and the specific heat ratio 7 = |. The speed of light is c = 10 V,,. The 
reconnection rates do not change noticeably when the ratio of the speed of light to the 
Alfven speed is increased to c = 100 as was done in [47|. The initial number densities 
are given by, 

ne = ni = no(l+sech^(^)]. (5.7) 



.5 VA , 

The electron and ion temperatures differ slightly, but are constant throughout the do- 
main, this gives the following electron pressure, Pg, 

Bg? (5.8) 



12^(0 '"0 
and ion pressure P, 

P, = -^b2^. (5.9) 
The electron and ion pressure balance the magnetic field which is given by 



Bo f2n\ . f2nx\ ny 
— — sm — — cos — ^ 

10 V i^;. y \u j V ^ 



B.v = =7^(;^)sin(^)cos(^) (5.11) 
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The magnetic field is in equilibrium with with the electron current Jz 



he = ^sech'('). (5.12) 



A VA 

The simulations are run at two resolutions, 128 x 64 and 512 x 256 using a second or third 
order discontinuous Galerkin method with third order TVD Runge Kutta time stepping. 
At a resolution of 512 x 256, 80 thousand time steps are taken and at a resolution of 128 x 



64, 20 thousand time steps are taken. In figure 14 low resolution solutions to the GEM 
challenge problem are plotted at f = 25/cOci- The second order solution shows stable 
island formation while in the third order method no island forms. In this simulation the 
TVB limiter constant M = is used to eliminate the formation of an unstable island in 
the third order method. Solutions to the GEM challenge problem are highly susceptible 
to bifurcation | [48| which is the formation of magnetic islands at or near the x-point (the 
center of the domain in this case). The development of these islands may be due to 
excessive dissipation applied at the x-point; however, the development of islands can be 
unpredictable and in the case of the 3rd order method, increased dissipation can actually 
eliminate the island. 

The islands are unstable and any small perturbation in the location or fields in the 
island will cause the island to slip and merge with one of the lobes. Small perturbations 
can arise from machine precision error and this is particularly true at locations near equi- 
librium where two nearly equal but opposing forces are added, the number that remains 
may depend significantly on the precision of the numbers used. It has been observed that 
increasing machine precision can change the direction that an unstable island slips off to. 
The second order solution shows the formation of a stable island presumably due to the 
extra dissipation of the second order method. 



In figure 15 the high resolution solutions are plotted at t=25/ cOd with the TVB limiter 
constant M = 25. At this resolution the second and third order methods produce very 
similar results. In this case no island formation is visible in either the second or third 
order schemes, however at time t = 30/ ood a very small unstable island forms in the 



3rd order method and combines with the left lobe before t = 35/ cOci- In figure 16 the 
reconnected magnetic flux ^ j\Bxj\dx along the x axis is plotted for several solutions. At 
grid resolution of 512 x 256 the second and third order methods produce similar results. 
At a resolution of 128 x 64 the third order method produces results which are in much 
better agreement with the high resolution results than does the second order method. 



In figure 17 the reconnected flux up to time i =60/ cOd is shown illustrating the satura- 
tion of the numerical solution beyond the published time interval in the GEM challenge 
results 1^. 



Cl- 



Flows appear late in time which are turbulent as shown in figure 18 at time ^=40/0;, 
The 3rd order method shows considerably more asymmetry than the 2nd order method, 
but asymmetry has begun to develop in the 2nd order solution. This asymmetry develops 
out of machine precision errors that are initiated in the current layer at the x-point from 
the balancing of source terms and fluxes. The only dissipation present is numerical so 
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Figure 14: GEM challenge comparing I/2I using the 2nd and 3rd order method on a 128x64 grid at time 
1 = 75 1 w^{. At this resolution an island forms in the 2nd order method and grows as the simulation progresses. 
Both methods use 3rd order TVD Runge-Kutta time stepping. 
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Figure 15: GEM challenge comparing \}z\ using the 2nd and 3rd order method on a 512x256 grid at time 
t = 25/ Wi-i. At this resolution the 2nd and 3rd order method look similar. Both methods use 3rd order TVD 
Runge-Kutta time stepping. Instabilities develop out of the fluid jets that emerge from the x-point when they 
collide with the slow moving fluids in the lobes. 
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Figure 16: Plot of reconnected flux vs time for 2nd and 3rd order spatial discretizations at resolutions of 128x64 
and 512x256. The reconnected flux differs substantially for the two methods at 128x64, though the 3rd order 
method better matches the high order solutions. At 512x128 the 2nd and 3rd order methods are in close 
agreement. 
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Figure 17: Plot of reconnected magnetic flux to time t = 60/cOci for a 512X256 solution. This figure shows that 
the two-fluid solution does saturate as would be expected with conducting walls and finite flux in the domain. 
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the second order method tends smooth out the errors produced by the finite precision 
more than the 3rd order method. As a result of the low dissipation in the 3rd order 
method these errors result in the excitation of unstable modes which eventually effect 
the macroscopic solution. The magnitude of electron momentum can differ by as much 
as 10 percent by moving from 64 bit numbers to 80 bit numbers at i =40/ cOd', fortunately 
these regions tend to be localized. 

Notice the pair of shocks in figure 18 the positions differ in the two solutions. The 
positions differ because the onset of the fast growth stage is slightly different for the two 
methods. As the fast growth begins, shocks form in the ion fluid as it is accelerated along 
the X axis. Eventually the two jets of shocked ion fluid collided and two shock that spans 
the y axis are formed which continues to propagate through the domain as the solution 
evolves. If the onset of fast growth differs slightly initially, then these shocks will appear 
at different locations for a given time. 




Figure 18: Plot of total ion momentum at t = 40/cOi.i at a resolution of 512 x 256 using the second and third 
order DG methods. 
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6 Discussion 

A discontinuous Galerkin method for the ideal 5-moment two-fluid plasma system is 
developed. A scalar model problem of the ideal two-fluid system is derived to illustrate 
the character of the full system. An analytic two-fluid solution to an electron acoustic 
square pulse in the linear regime is derived and the numerical solution using the fully 
non-linear two-fluid system is calculated showing convergence for the 2nd and 3rd order 
discontinuous Galerkin methods. The algorithm is benchmarked against the two-fluid 
electromagnetic shock originally published in The 2nd and 3rd order algorithms 
are tested on the GEM challenge magnetic reconnection problem and produces results 
comparable to those generated by particle codes, hybrid codes, and a Hall MHD code in 
p9| . The discontinuous Galerkin method offers a straight forward method for solving the 
ideal 5-moment two-fluid system. This same algorithm could be applied to 10 and higher 
moment two-fluid systems. The algorithm can be easily generalized to three dimensions 
and general geometries. 
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